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ABSTRACT 


The  area  of  the  Bering  Sea  and  Alaska  has  been  studied 
in  terms  of  shear-velocity,  density  and  compress ional-velo- 
city  structure  by  applying  a  generalized  linear  inversion 
method  to  fundamental-mode  Rayle igh-wave  group-velocity 
dispersion  relationships  in  the  period  range  from  10  to 
100  sec. 

Group  velocity  dispersion  relationships  in  the  area 
have  been  obtained  by  applying  the  phase-matched  filtering 
technique  (-He-rrin  and  Goforth,  1977)  to  digitally  recorded 
surface-wave  data.  Corrections  for  instrument  response  and 
the  sphericity  of  the  Earth  were  applied  to  the  dispersion 
observations.  A  new  exact  analytical  method  for  the  compu¬ 
tation  of  Rayle igh-wave  phase-velocity  partial  derivatives 
with  respect  to  Earth  parameters  has  been  formulated.  With 
the  phase-velocity  partial  derivatives  determined,  the  group 
velocity  partial  derivatives  were  computed  by  use  of  the 
fast  and  accurate  method  of  Rod i  eb-  a 1 .  (1975),  and  were 

successfully  incorporated  into  a  generalized  linear  inversion 
method . 


IV 


The  study  area  has  been  found  to  consist  of  three 
physiographic  provinces  and  the  structure  of  the  three 
regions  has  been  estimated  as  follows:  In  continental 

i  ■.  - 

Alaska,  the  crustal  thickness  is  43  t  3  km,  and  a  low 
velocity  zone  extends  from  a  depth  of  about  113  km  to  about 
213  km.  In  the  Bering  Shelf  region,  the  depth  to  the 

t  ct 

Mohorovicic  discontinuity  is  28  4  km,  and  a  low  velocity 

zone  ranges  in  depth  from  about  108  km  co  about  213  km.  In 

-i  'i  •••' 

the  Aleutian  Basin,  the  thickness  of  the  crust  is  18  -  4  km, 
and  a  low  velocity  zone  extends  from  a  depth  of  about  60  km 


to  about  220  km. 
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INTRODUCTION 


Seismic  surface  waves  have  drawn  much  of  seismologists' 
attention  in  the  last  few  decades,  partly  due  to  their  ability 
to  sample  the  outer  part  of  the  earth  over  long  paths  and  to 
provide  information  on  average  structures  for  parts  of  the 
earth  which  are  not  readily  accessible  to  body  wave  studies; 
and  partly  due  to  their  usefulness  in  discrimination  between 
natural  and  artificial  seismic  sources,  and  in  estimation  of 
the  source  parameters . 

The  objective  of  this  study  is  to  further  our  understand¬ 
ing  of  the  structure  of  the  earth  by  investigating  the  crustal 
and  uppermost  mantle  structure  of  the  Bering  Sea  and  of 
Alaska.  In  order  to  accomplish  this  goal,  group  velocity 
relationships  have  been  determined  by  applying  the  phase- 
matched  filtering  technique  (Herrin  and  Goforth  1977)  to 
observed  surface  wave  data.  In  the  course  of  inverting  the 
group  velocity  dispersion  relationships,  in  order  to  obtain 
the  shear-wave  velocity,  density  and  compressional-wave 
velocity  structure  of  the  area,  a  new  method  for  the  computa¬ 
tion  of  Rayleigh-wave  phase-velocity  partial  derivatives  with 
respect  to  earth  parameters  was  formulated.  With  the  phase- 
velocity  partial  derivatives  determined,  the  fast  and 
accurate  method  of  Rodi  et  al_.  (1975)  was  used  for  computing 
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the  group-velocity  partial  derivatives.  In  conjunction  with 
this  new  method,  a  generalized  linear  inversion  technique 
was  used  for  the  inversion  of  the  group  velocity  dispersion 
relationships . 

The  area  of  the  Bering  Sea  and  Alaska  was  chosen  to 
be  investigated  for  the  following  reasons:  (i)  This  area 
bears  geophysical  and  tectonic  importance  in  that  it  is  loca¬ 
ted  along  the  plate  boundary  between  the  Pacific  and  North 
American  plates  and  is  adjacent  to  the  boundary  between  the 
North  American  and  Eurasian  plates.  (ii)  This  area  is  a 
region  of  primary  significance  for  seismic  surveillance  from 
the  United  States'  point  of  view  (Evernden,  1969).  (iii)  This 
area  may  be  unusual  in  that  the  crust  of  the  Aleutian  Basin 
region  may  be  undergoing  the  process  of  the  continentaliza- 
tion  of  an  oceanic  crust  (Shor,  1964) .  (iv)  Finally,  this 
area  has  not,  to  this  author's  knowledge,  been  studied 
before  in  detail  and  on  a  regional  basis . 

It  was  decided  to  work  on  group  velocities  rather  than 
phase  velocities  with  the  following  reasons:  (i)  As  men¬ 
tioned  above,  the  phase-matched  filtering  technique,  which  is 
an  accurate  method  for  determining  group-velocity  dispersion, 
was  readily  available.  (ii)  It  is  well  known  that  group  and 
phase  velocity  observations  provide  basically  the  same  infor¬ 
mation  about  the  earth  model  (Der,  et_  al . .  1970;  Wiggins, 
1972) .  Pilant  and  Knopoff  (1970)  have  shown  that  inver¬ 

sion  of  phase  velocities  eliminates  one  degree  of 


non-uniqueness  compared  wich  inversion  of  group  velocities 
since  group  velocities  are  a  derivative  property  of  phase 
velocities.  However,  as  pointed  out  by  Yu  and  Mitchell 
(1979) ,  due  to  the  very  nature  of  the  group  velocity  as  a 
differential  of  the  phase  velocity,  small  perturbations  in 
phase  velocity  show  up  as  larger  variations  in  group  velo¬ 
city.  Thus  the  group  velocity  dispersion,  when  inverted, 
should  produce  a  more  detailed  earth  model.  Braile  and 
Keller  (1975)  have  stated  that  the  inversion  of  group  velo¬ 
cities  did  not  appear  to  have  degraded  the  results  in 
comparison  with  those  from  the  inversion  of  phase  velocities, 
(iii)  The  aperture  of  the  Alaskan  Long  Period  Array  (ALPA) 
used  as  the  recording  station  may  be  too  small  to  yield 
reliable  phase  velocity  information. 

Preliminary  inspection  of  data  showed  that  the  quality 
of  Love-wave  data  was  inferior  to  that  of  Rayleigh-wave  data. 
It  was  also  found  that  Love-wave  dispersion  contributes  less 
information  about  earth  structure  than  Rayleigh-wave  disper¬ 
sion  (Der  et  al . ,  1970;  Braile  and  Keller,  1975).  Since,  in 
general,  the  maxima  of  the  Rayleigh-wave  partial-derivative 
curves  are  larger  and  narrower  than  those  of  the  corresponding 
Love-wave  curves,  Rayleigh  waves  give  better  resolution  than 
Love  waves.  For  the  above  reasons,  it  was  decided  to  use 
only  Rayleigh  waves  in  this  study. 

Inspection  of  results  of  the  multiple  filter  analysis 
(Dziewonski  et  al.,  1969),  showed  that  there  was  not  much 


energy  present  in  higher  modes.  Under  this  circumstance, 
efforts  to  use  higher  modes  may  introduce  erroneous  rather 
than  useful  information.  Therefore  it  was  decided  to  work 
with  only  fundamental -mode  Rayleigh  waves. 


DATA 


The  data  used  in  this  study  were  collected  at  the 
Alaskan  Long  Period  Array  (ALPA) .  ALPA  was  selected  as  a 
recording  station  for  the  following  reasons:  (i)  This  seis¬ 
mic  array,  which  is  located  just  north  of  Fairbanks,  Alaska, 
was  probably  the  best,  reliably  equipped  set  of  seismographic 
stations  around  the  study  area  and  data  from  it  are  suitable 
for  accomplishing  the  objective  of  this  research  in  that  the 
array  is  located  at  one  end  of  the  region  of  interest,  while 
the  locations  of  the  seismic  sources  range  from  near  the 
array  to  the  Komandorsky  Islands.  (ii)  The  quality  of  the 
data  recorded  at  this  array  is  known  to  be  good.  It  is  an 
array  with  an  aperture  of  about  80  km  and  was  installed  with 
the  knowledge  and  experience  gained  from  the  operation  of  the 
large  aperture  seismic  array  near  Billings,  Montana  (LASA) . 
(iii)  Since  the  data  is  digitally  recorded  on  magnetic  tape, 
there  is  no  possibility  of  error  incurred  from  the  digitiza¬ 
tion  of  analog  data.  (iv)  Finally,  the  data  were  available 
for  use  in  this  research. 

The  location  of  the  recording  stations  and  the  epicenters 
of  the  earthquakes  used  in  this  study  are  shown  in  Fig.  1.  The 
relative  amplitude  and  phase  responses  of  the  seismograph  sys¬ 
tem  are  shown  in  Figs.  2  and  3,  respectively.  The  parameters 
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Figure  1:  Map  of  the  Study  Area  with  Locations  of  the 
graphic  Station  (ALPA)  and  Epicenters  of  the 
Seismic  Events  Studied  Shown.  Mercater  Prol 


of  the  earthquakes  used  are  listed  in  Table  1.  The  informa¬ 
tion  on  source  parameters  is  from  the  Preliminary  Determination 
of  Epicenters  (PDE)  published  by  either  the  National  Oceanic 
and  Atmospheric  Administration  or  the  United  States  Geological 
Survey  National  Earthquake  Information  Service .  The  vertical- 
component  seismograms  of  the  events  studied  are  shown  in 
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METHOD 


1.  Phase-matched  Filtering 

For  the  determination  of  the  dispersion  relationships 
of  surface  waves,  Dziewonski  et  al.  (1969)  have  developed  a 
method  called  multiple  filter  analysis.  Realizing  the  lack 
of  resolution  in  the  time  domain  of  the  above  method,  a  more 
capable  method  called  phase-matched  filtering  has  been  devel¬ 
oped  and  proved  effective  in  more  recent  studies  by  Herrin 
and  Goforth  (1977) .  The  phase-matched  filtering  technique 
was  used  in  this  study  in  order  to  separate  the  effects  of  the 
multipa thing  from  those  of  primary  arrivals.  In  using  the 
phase-matched  filter,  a  dispersion  curve  obtained  from  the 
multiple  filter  analysis  was  used  as  the  initial  input  to  the 
iterative  procedure  of  the  phase-matched  filtering  technique. 

Since  the  original  papers  give  full  accounts  of  the 
above  methods,  only  a  brief  description  of  the  methods  is 
given  below.  The  multiple  filter  analysis  technique  displays, 
as  its  product,  the  instantaneous  amplitude  of  a  seismogram 
as  a  function  of  both  period  and  group  velocity  (time) .  When 
the  instantaneous  amplitudes  are  contoured,  a  line  connecting 
the  maximum  amplitudes  for  specified  periods,  which  is  rep¬ 
resented  by  the  ridge  crest  of  the  contour  diagram,  gives 
the  group  velocity  dispersion  relationships  for  the  analyzed 
seismogram. 
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In  the  phase-matched  filtering  technique,  the  Fourier 
phase  of  the  filter  is  made  equal  to  that  of  a  given  seismo¬ 
gram  with  the  possibly  existing  multipath  effects  removed. 

The  group  delays  for  the  initial  input  dispersion  curve  are 
computed.  The  phase  of  the  filter  is  calculated  by  inte¬ 
grating  the  above-obtained  group  delays.  By  crosscorrelating 
the  filter  and  the  given  seismogram  and  windowing  the  cross¬ 
correlation  function,  the  phase  of  the  filter  is  made  equal 
to  that  of  the  seismogram  after  a  certain  number  of  iterations. 
The  group  delays  of  the  filter  are  corrected  as  the  phase  of 
the  filter  is  modified  in  each  iteration.  The  group  velocity 
is  then  computed  from  the  final  group  delays  of  the  phase - 
matched  filter. 

2.  Division  of  the  Study  Area  into  Different  Provinces 

The  study  area  was  divided  into  three  physiographic 
provinces--continental  Alaska,  the  Bering  Shelf,  and  the 
Aleutian  Basin.  This  classification  was  chiefly  based  on 
topographic  and  bathymetric  data  of  the  area  and  is  shown 
in  Fig.  1.  The  boundary  between  the  Bering  Shelf  and  the 
Aleutian  Basin  was  drawn  along  the  1000-fathom  contour 
line  of  bathymetric  data.  One  purpose  of  this  study  was 
to  study  the  structure  of  the  three  provinces  and  to  deter¬ 
mine  the  boundary  between  them  if  they  are  different  in 
structure . 

Although  it  appeared  from  the  topographic  and  bathymetric 
data  that  continental  Alaska  and  the  Bering  Shelf  are 


continental  and  the  Aleutian  Basin  is  oceanic  in  crustal  and 
uppermost  mantle  structure,  there  was,  as  pointed  out  in  the 
introduction,  a  possibility  of  the  Aleutian  Basin's  being 
continental  in  nature.  As  an  attempt  to  resolve  this 
problem,  the  following  method  was  used. 

By  inverting  the  group  velocity  dispersion  relationship 
obtained  from  the  event  0509,  whose  epicenter  is  located  with¬ 
in  continental  Alaska,  a  structure  for  continental  Alaska  was 
determined.  It  was  recognized  that  the  great  circle  path  for 
the  event  1410,  which  is  located  near  the  western  edge  of  the 
Bering  Shelf,  crosses  the  boundary  between  continental  Alaska 
and  the  Bering  Shelf  at  a  nearly  right  angle .  The  proportion 
of  the  lengths  of  continent  and  shelf  along  the  great  circle 
path  was  measured  on  a  globe.  By  subtracting  the  group  velo¬ 
city  dispersion  proportional  to  the  continental  path  from  the 
observed  group  velocity  dispersion  for  the  event  1410,  the  ob¬ 
served  group  velocity  dispersion  relationship  for  the  pure  shelf 
path  was  derived.  By  inverting  the  resultant  group  velocity 
dispersion,  the  structure  of  the  shelf  region  was  determined. 

Similarly,  it  was  found  that  the  great  circle  path  for 
the  event  2224,  which  is  located  near  the  western  end  of  the 
Aleutian  Basin,  crosses  the  boundary  between  continental  Alaska 
and  the  Bering  Shelf  and  the  boundary  between  the  latter  and 
the  Aleutian  Basin,  both  nearly  perpendicularly.  The  group 
velocity  dispersion  relationships  determined  above  for  con¬ 
tinental  Alaska  and  the  Bering  Shelf  were  subtracted  in 
proportion  to  their  segment  lengths  of  the  great  circle  path 


from  the  observed  group  velocity  dispersion  of  event  2224 
to  derive  the  observed  group  velocity  dispersion  relationship 
for  the  pure  Aleutian  Basin  path. 

Before  committing  much  effort  to  find  the  detailed 
structure  of  the  study  area,  reasonableness  of  the  above  line 
of  thinking  was  tested  by  applying  the  above  procedure  with 
average  earth  models  of  Dziewonski  et  al.  (1975)  to  the  ob¬ 
served  data  corrected  for  instrument  response.  The  results 
were  encouraging,  and  the  above  division  of  the  study  area  was 
adopted  for  subsequent  studies.  The  proportion  of  lengths  of 
different  provinces  along  the  great  circle  path  for  the  events 
used  in  this  study  is  shown  in  Table  2. 

3.  Corrections  to  the  Observed  Data 

Two  types  of  corrections  were  applied  to  the  group  velo¬ 
city  dispersion  relationships  determined  from  the  phase-matched 
filtering  technique.  One  is  for  the  instrument  phase  res¬ 
ponses  and  the  other  is  for  the  sphericity  and  gravity  of  the 
earth . 

For  the  instrument  correction,  the  following  relationship 
was  used: 

'■> 

where  and  are,  respectively,  the  group  delay  and 

phase  response  of  the  instruments  for  an  angular  frequency  o. 
The  phase  responses  of  the  instruments  as  shown  in  Fig.  3 
were  used  for  this  purpose. 
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table  2 

PROPORTION  OF  PATH  LENGTHS  IN  DIFFERENT  PROVINCES 


Event 

Name 

Proportion  of  Distances 
Continental  Shelf 

Oceanic 

1410 

0.50  0.50 

0.00 

1519 

0.51  0.30 

0.19 

2005 

0.39  0.27 

0.34 

2224 

0.25  0.40 

0.35 

During  the  summer  of  1972,  the  seismograph  filter- 
amplifiers  were  replaced  with  new  units  having  different 
characteristics .  Two  different  response  curves  for  the  cor¬ 
responding  filter-amplifier  systems  are  shown  in  Fig.  3. 

Events  0509  and  2005  were  recorded  with  the  old  instruments 
while  the  rest  of  the  events  were  recorded  with  the  new 
instruments . 

The  observed  group  velocities  as  obtained  from  the  phase- 
matched  filtering  technique  represent  the  dispersion  of  the 
gravitating  spherical  earth,  while  the  calculation  of  the  group 
velocities  in  the  process  of  inversion  is,  for  computational 
convenience,  made  for  a  nongravitating  plane-layered  earth 
model.  Therefore  it  is  necessary  to  convert  the  observed 
group  velocity  dispersion  to  one  corresponding  to  the 
nongravitating  plane-layered  earth. 

For  Love  waves  it  is  possible  to  apply  a  suitable  trans¬ 
formation  to  the  plane-layered  earth  model  to  simulate  the 
effects  of  sphericity  (Gerver  and  Kazhdan,  1968;  Biswas  and 
Knopoff,  1970).  For  Rayleigh  waves,  however,  a  similar  trans¬ 
formation  is  difficult  to  achieve  and  an  empirical  correction 
must  be  applied.  Such  a  correction  has  been  derived  by  Bolt 
and  Dorman  (1961)  and  has  been  widely  used.  After  fifteen 
years,  North  and  Dziewonski  (1976),  based  on  studies  with 
additional  and  presumably  better  earth  models,  have  improved 
the  formula  of  Bolt  and  Dorman  (1961) . 

Since  the  effects  of  sphericity  and  gravity  upon  group 
velocity  are  much  less  than  those  upon  phase  velocity,  as 


pointed  out  by  Bolt  and  Dorman  (1961)  and  North  and 
Dziewonski  (1976) ,  a  quarter  of  the  difference  between  the 
phase  velocities  of  gravitating  spherical  earth  and  nongravi¬ 
tating  plane-layered  earth,  as  shown  in  Table  1  of  North  and 
Dziewonski  (1976)  was  used  in  the  present  study  as  a  somewhat 
arbitrary  but  reasonable  difference  between  the  group  velo¬ 
cities  of  the  above  two  representations  of  the  earth  (see 
Figs.  4  to  6  of  Bolt  and  Dorman,  1961).  For  the  regions  of 
continental  Alaska,  the  Bering  Shelf  and  the  Aleutian  Basin, 
the  velocity  differences  of  the  earth  models  ANDES,  PLATFORM 
and  PEMOCD,  respectively,  of  North  and  Dziewonski  (1976)  were 
used.  Since  it  was  found  that  the  velocity  differences  of 
North  and  Dziewonski  are  smoothly  varying  quantities  with 
period,  their  values  were  interpolated  or  extrapolated  as 
necessary  in  the  present  study. 

4.  Uncertainty  in  the  Observed  Data 

Der  et  al.  (1970)  used  the  uncertainty  principle  of 
Bendat  (1958)  in  assessing  the  uncertainties  in  their  group 
velocity  data  obtained  by  the  multiple  filter  analysis  tech¬ 
nique.  Judging  that  the  multiple  filter  analysis  gives  an 
accuracy  greater  than  that  from  the  uncertainty  principle, 
they  took  one-fifth  of  the  error  estimated  from  the  uncer¬ 
tainty  principle  as  their  error  of  measurements.  They  gave 
0.01  km/sec  as  a  lower  error  limit;  0.03  km/sec  was  adopted 
as  the  uncertainty  of  the  group-velocity  data  by  Braile  and 
Keller  (1975)  without  an  explicit  explanation,  and  by  Yu  and 
Mitchell  (1979)  who  used  the  multiple  filter  analysis  method. 
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Although,  in  the  present  study,  the  phase-matched  fil¬ 
tering,  which  is  considered  to  be  much  more  accurate  than  the 
multiple  filter  analysis,  is  used,  0.03  km/sec  was  taken  as 
the  maximum  error  in  the  group-velocity  data.  This  value  is 
adopted  in  order  to  accomodate  the  possible  errors  in  applying 
the  instrument  and  sphericity-gravity  corrections. 

5.  Method  for  the  Computation  of  Rayleigh-wave  Phase-velocity 
Partial  Derivatives  with  Respect  to  Earth  Parameters 

In  the  inversion  of  group  velocity  dispersion  observations 
it  is  necessary  to  compute  the  partial  derivatives  of  group 
velocities  with  respect  to  various  earth  parameters.  Since  a 
fast  and  accurate  method  has  been  developed  for  computing  group 
velocity  partial  derivatives  when  corresponding  phase  velocity 
partial  derivatives  are  known  (Rodi  et:  al. ,  1975),  an  effort 
was  made  in  the  present  study  to  formulate  a  simple  new  method 
for  the  computation  of  exact  phase-velocity  partial  deriva¬ 
tives.  The  following  paragraphs  are  the  outcome  of  this 
effort . 

In  an  earth  model  in  which  the  earth  consists  of  many 
elastic,  homogeneous  and  isotropic  parallel  layers  on  a  homo¬ 
geneous  and  isotropic  half-space,  the  earth  properties  may  be 
defined  by  longitudinal  and  transverse  wave  velocities,  den¬ 
sity,  and  thickness  of  the  component  layers.  Various  methods 
have  been  employed  in  computing  the  phase-velocity  partial 
derivatives  with  respect  to  the  above  model  parameters. 

Dorman  and  Ewing  (1962),  followed  by  Brune  and  Dorman 
(1963) ,  calculated  the  changes  in  phase  velocity  due  to  the 


perturbation  of  each  physical  parameter  while  retaining  the 
remaining  parameters  constant.  This  method  is  very  time- 
consuming.  On  the  basis  of  Jeffreys'  (1961)  suggestion  that 
Rayleigh's  principle  can  be  used  to  find  expressions  for  the 
effects  of  small  changes  of  the  elastic  properties  on  the 
phase  velocity,  Anderson  (1964)  and  Takeuchi  et  al.  (1964), 
followed  by  Harkrider  (1968)  and  Anderson  and  Harkrider  (1968) , 
approached  the  problem  by  using  the  energy  integral  technique. 
This  technique  is  quite  complex.  McEvilly  (1964)  has  used  a 
combination  of  the  above  two  methods.  Bloch  ej:  al.  (1969), 

Der  et  al_  (1970),  Der  and  Landisman  (1972)  and  Knopoff 
(1972)  have  computed  phase-velocity  partial  derivatives  but 
they  have  not  described  the  algorithm  they  employed.  Novotny 
(1970)  derived  exact  expressions  for  Love-wave  phase  velocity 
partials  in  a  different  way  by  taking  advantage  of  Thomson- 
Haskell  matrices  (Thomson,  1950;  Haskell,  1953). 

In  the  rest  of  this  section,  following  Novotny  in  Love 
waves,  a  new  method  of  computing  the  Rayleigh-wave  phase- 
velocity  partial  derivatives  with  respect  to  the  parameters 
of  the  medium  is  presented.  Convenient  formulas  for  the 
Haskell  layer-matrix  derivatives,  which  are  necessary  in  the 
computation  of  the  phase-velocity  partial  derivatives,  are 
given  in  the  Appendix.  Employing  double  precision  in  computer 
programming,  the  use  of  the  Haskell  layer-matrix  method  has 
not  shown  any  numerical  difficulties  (Thrower,  1965;  Dunkin, 
1965;  Gilbert  and  Backus,  1966,  among  others)  in  the  period 
range  of  10  to  100  sec,  as  used  in  this  study. 


Harkrider  (1964)  has  derived  a  matrix  equation  for  the 
dispersion  of  Rayleigh  waves,  which  may  be  expressed  as  follows 


F (C,  o,  dm,  /S»,  A.,  A.)  =  N K  i  L*rt*  -  T’Cg’n  -  L* H)  U) 
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Here  c*?  is  the  angular  frequency,  c  phase  velocity,  oC  and 
|3  velocities  of  compressional  and  shear  waves ,  /  density,  d 
thickness  of  a  layer; 

Xy^  -  2(§m./c)  , 

=  4l  dm.  vim  ,  Qm  -  dl  dm  dip  m 


where  d.  represents  the  wave  number; 
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where  is(-i),/x; 

^  =  Q-n-i  '  ■  ' 

with  designating  the  Haskell's  layer  matrix  for  the  mth 
layer;  and  subscripts  m,  n  and  o  refer  to  the  mth  solid  layer, 
half-space  and  liquid  surface  layer,  respectively.  In  order 
to  minimize  typographical  requirements,  a  compound  subscript, 
e.g.  oU n  ,  is  used  where  a  subscripted  subscript  is  appropriate. 
An  asterisk  superscript  is  used  to  transform  an  imaginary 
quantity  into  a  real  one  such  that  for  an  arbitrary 

imaginary  quantity  £  . 

Since  T*~°  for  the  earth  model  in  which  the  liquid 
surface  layer  does  not  exist,  considering  (2)  is  sufficient 
for  earth  models  both  with  and  without  ocean.  From  (2), 
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where  X  designates  one  of  the  earth  parameters  oC ,  ^  >  j*  and  (£• 
If  we  define  two  product  matrices  and  such  that 
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For  the  particular  cases  where  ?»i®/  and  7n*n-)  ,  respectively, 
the  second  equation  of  (6)  becomes 

ov  l  1  9Ai  Jij 


3  r  9 


"Mix  M-- 

>1-1  *"  9  «n-l  -*  t  J 


It  is  noted  from  (2)  and  (3)  that  in  order  to  compute 
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It  is  also  noted 


,  in  particular,  we 


and  terms  necessary  to  compute 


Since 


and  '7 —  are  simple  to  calculate  they  are  not  discussed  here. 

From  the  above  discussion,  it  becomes  clear  that  in 

order  to  compute  phase-velocity  partial  derivaties  with  res- 

j  n  _  902»t),f  3 

pect  to  model  parameters,  we  need  only  to  compute  — L,  — - — 

(?Xm  9  " 


and 


The  expressions  for  these  terms  are  given  in 
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the  Appendix. 

The  merit  of  the  present  method  may  lie  in  the  simplicity 
of  the  algorithm.  Equations  (4)  and  (6)  in  connection  with 
dispersion  equation  (2)  explain  essentially  all  procedures 
required.  Partial  derivatives  of  layer-matrix  elements  given 
in  the  Appendix  provide  expressions  needed  in  the  computation 
of  (4)  and  (6) . 

6.  Generalized  Linear  Inversion 

The  generalized  linear  inversion  method,  which  may  be 
viewed  as  a  discrete  specialization  of  the  general  formulation 
of  the  Backus-Gilbert  inversion  technique  (Backus  and  Gilbert, 
1967,  1968,  1970),  has  been  applied  to  geophysical  problems 
in  many  previous  studies  (Smith  and  Franklin,  1969;  Der  et  al .  , 
1970;  Parker,  1970;  Wiggins,  1972;  Jackson,  1972;  Ward  et  al . , 
1973;  Braile  et  al . ,  1974;  Braile  and  Keller,  1975;  Pedersen, 
1977,  among  others).  The  inversion  scheme  used  in  this  study 
closely  follows  the  algorithm  described  by  Wiggins  (1972)  and 
Jackson  (1972).  However,  an  outline  of  the  method  along  with 
some  points  relevant  to  the  present  study  are  given  below. 

In  general  geophysical  inversion  problems,  we  have  a 
nonlinear  system  dealing  with  n  observations  and  m  model  para¬ 
meters.  In  the  following  we  shall  use  the  subscript  i  to 
represent  the  ith  observation  and  the  subscript  j  to  designate 
the  j th  parameter.  We  shall  use  a  small  letter  to  represent 
a  scalar,  a  small  letter  with  a  bar  underneath  it  a  vector, 
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and  a  capital  letter  a  matrix. 

The  earth  model  we  assume  must  be  capable  of  providing 
us  with  a  functional  relationship  between  the  model  parameters 
and  the  calculated  values.  The  model  parameters  will  be 
represented  by  a  vector  such  that 

X  =  (oC,  /3,  P.  d) 

where  ,  (3  ,  P  and  d  are  functions  of  depth.  The  calcu¬ 

lated  values  represent  the  group  velocity  and  will  be  designated 
by  .  We  expand  the  problem  functional  "*$(&)  in  a 

Taylor  series  about  an  initial  model  Xe  : 


?lrr.  V  * 

— u-u) 


- % (M 


Wits) 

5Xj  /*.* 


■  (!-*,)  +  -kiZ  (1) 

KsX. 


J> 

where  is  the  functional  in  the  second  or 

higher  (■£  M)  order  terms  in  -  £0)  .  We  linearize 

the  problem  by  retaining  the  linear  part  and  neglecting  the 
nonlinear  part  of  (7) .  By  setting 

SXj  =  Xj  -  Xoj 

a-  -  Hk 
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(7)  becomes 

£v  -  A  £X  m 

Hence  we  have  a  linear  system  in  n  equations  and  m  unknowns 
and  the  systems  Matrix  A  is ,  in  general,  not  square. 

Equation  (8)  may  be  viewed  as  follows:  Xj  represents 
the  true  earth,  X0j  the  initial  model,  <fX j  the  correc¬ 
tion  to  be  made  to  the  initial  model;  (X)  represents 

the  observed  group  velocity  at  the  ith  period, 
the  group  velocity  for  the  ith  period  computed  from  the  ini¬ 
tial  model,  6H.  the  difference  in  group  velocity  between 
the  model  and  the  true  earth.  It  is  noted  that  the  expansion 
of  the  problem  functional  in  Taylor  series  is  based  on  the 
assumption  that  the  group  velocity  is  a  smoothly  varying 
function  of  the  parameters.  Also  the  linearization  is  justi¬ 
fied  only  for  small  values  of  the  parameter  corrections  <fXj  . 
If  SXj  must  be  large  to  satisfy  the  observations,  then  the 
results  must  be  checked  by  expanding  the  functional  about  a 
new  initial  model. 

For  simplicity  in  notation,  we  replace  for  <fX  and 
\  for  £>  ,  while  keeping  in  mind  that  X  and  actually 
represent  <£X  and  in  ,  respectively.  Then  (9)  becomes 

A  *  =  H  do) 

If  we  can  find  A  ,  the  exact  inverse  of  A  ,  then 
we  can  solve  the  system  for  X  .  But  since  A  is ,  in  general, 
not  square,  and  possibly  singular  even  when  it  is  square,  the 
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inverse  of  A  does  not  exist.  The  generalized  linear  inver¬ 
sion  is  a  way  of  avoiding  this  difficulty  by  finding  the  best 
possible  approximation  to  a  matrix  of  quantities,  which  may 
be  conceivably  called  the  exact  inverse  of  A  and  designated 
by  A*  .  This  procedure  may  be  conceptualized  in  the 
following  way. 

First  consider  the  adjoint  system  of  (10) : 

ATX'*V  (W 


where  the  superscript  T  represents  the  transpose  of  the 
matrix.  The  combination  of  (10)  and  (11)  into  an  identity 
gives 


[ 0  ;  A ' 

■  X'- 

II 

- - — 
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o 

X  . 

1 

3* 

(12) 


It  is  noted  that  the  matrix  in  (12) ,  whose  blocks  consist  of 
A  ,  /\T  and  O  »  is  symmetric  and  its  size  is  (n+n)  x  (n+m). 

We  shall  call  this  symmetric  matrix  £  .  Let  us  consider 
the  eigenvalue  problem  for  $  : 

$  w  =A  w  (13) 

The  component  of  W  may  be  separated  into  contributions  from 
the  71 -dimensional  data  space  of  group  velocities  and  the 
TJt -dimensional  parameter  space  such  that 

WT=  (U,,  Uit  ■  ■  ■ ,Un .  V, .  •  •  •  Vm)  (14) 


Then  (13)  may  be  expressed  as 
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V 
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That  is , 

A  y  =  Au 

Atu  =  A  V 

From  (16)  we  obtain 

ata  t  =  A2  y 

AA’S 

Let 

/(  =  A*  ;  BSATA;  CMA7 

Then  (17)  becomes 
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(16) 


(17) 

(18) 


B  y  =  /*  y 

c  y  =/<  y 

From  either  of  the  standard  eigenvalue  problems  (19) ,  in 
connection  with  (16) ,  we  can  find  \  ,  U  and  V  . 

Since  0  and  C  are  symmetric,  Ui  and  Vj  form  sets 
of  orthogonal  vectors.  Let 

g  =  ^iW(6)  =  SLamA  CC)  (20) 

Then 


$  <  TfUrfL  (  72  ,  m) 


(21) 
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We  normalize  (J±  and  V^-  and  define  matrices  (J^  and  such 
that  the  columns  consist  of  >  and  Vj  ,  j-  i  •  $ 

respectively,  and  define  a  diagonal  matrix  Ag  whose  diagonal 
elements  are  ,  -i.  -  1  ,  £  arranged  in  decreasing  order. 

Then  we  realize  the  following: 

(i)  Always 

U?  -  J*  •'  V  V,  - 1,  «*> 

where  a  subscript  of  1  indicates  the  size  of  the  identity 
matrix.  Equation  (22)  is  true  since  the  multiplication  in  (22) 
always  involves  the  inner  product  of  orthonormal  eigenvectors . 

(ii)  When  ^  -  n  < 7H  ,  i.e.,  the  system  is  under-determined, 


aa  -Ir-l 


(23) 


and  this  condition  guarantees  the  existence  of  a  solution 
(Jackson,  1972) . 

(iii)  When  ^ *  m  <  n  ,  i.e.,  the  system  is  over-determined, 

<M) 

and  this  condition  guarantees  the  uniqueness  of  the  solution, 
if  a  solution  to  the  problem  exists  (Jackson,  1972) . 

(iv)  When  £  =  7t  =  Tn  ,  i.e.,  the  system  is  even-determined, 

u& u;  *  ^ 

and  this  condition  guarantees  the  existence  and  uniqueness  of 
a  solution. 

With  ur  v,  and  so  defined,  we  obtain  from  the 
first  equation  of  (16)  that 

A*U?A^V|T  (25) 


A. 


Substituting  (25)  into  (10)  gives 


U?A?V/x=* 

from  which  we  obtain 

H  H  u$  - 

Let 

R  =  V,  Vf 

Then 

R  *  =  Vg  Aj  Uj  H  (28) 

Let 

X  =  Rx  (i?) 

Then 

x  =  Vg  A*  Uf  a  (3°) 

From  (24)  and  (28),  it  is  noted  that  if  ^s»7n<n  ,  then  ft 
becomes  Iy*  and  (28)  gives  a  unique  solution  X  ,  which  is 
identical  to  X  ;  that  if  £+jh,  then  the  degree  of  likeness 
of  R  to  I  gives  the  degree  of  uniqueness  of  our  estimated 
solution  X  .  Due  to  this  property  of  R  ,  it  is  called  the 
resolution  matrix,  and  its  rows  are  called  resolving  kernels. 
Let 

Then  (30)  becomes 


(31) 
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from  which  it  is  easily  understood  that  J-J  bears  the  nature 
of  the  inverse  of  A  ,  through  which  we  could  find  the  best 


estimate  of  the  solution  £  .  Due  to  this  characteristic  of 

,  it  is  called  the  generalized  linear  inverse  of  A  • 

From  the  viewpoint  of  the  above  discussion,  (10)  may  be 

viewed  in  the  following  two  ways: 

•x 

n 

l«U? 

o 

(33) 

A  i  =5 

(34) 

where  the  superscripts  o  and  p  represent  the  observed  and 

predicted  (calculated  from  the  earth  model)  data,  respectively. 

From  this  point  of  view,  (32)  may  be  expressed  as 

2  •  H  f 

(35) 

Substituting  (35)  into  (34)  gives 

v  =  a  h  y 

(36) 

Again  substituting  (25)  and  (31)  into  (36)  gives 

s'- W 

Let 

(37) 

d  =  u,  v; 

(38) 

Then 

f  =  d  v 

(37) 
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From  (23)  and  (37) ,  it  is  noted  that  if  ^  =  71  <  791,  then  D 
becomes  In  and  (37)  tells  us  that  the  estimated  solution 
is  exact  in  the  sense  that  the  predicted  and  observed  data 
are  identical;  that  if  £  4  71 ,  then  the  degree  of  likeness  of 
D  to  I  gives  the  degree  of  independence  of  the  observed 
data.  Due  to  this  characteristic  of  P  ,  it  is  called  the 
information  density  matrix. 

7 .  Propagation  of  Errors 

Lastly  we  wish  to  know  how  the  uncertainty  of  our  obser¬ 
vations  affects  the  solution  of  the  problem  presented  by  (10) , 
which  is  an  error-free  system.  Since  we  have  a  certain  amount 
of  uncertainty  or  error,  £  ,  in  the  observed  data,  (10)  may 
be,  in  this  system  with  erroneous  data,  expressed  as 

§  =  A  x  -  %  (4o) 

-\ 

Then  the  best  estimate,  X  ,  of  the  true  solution  X  ,  is  the 
solution  which  minimizes  the  square  error 

€7§  =  (41) 

In  general,  different  data  may  have  different  reliability 
and/or  different  units.  In  order  to  take  this  into  account, 
we  wish  to  give  different  weights  to  different  data.  From 
the  fundamental  theory  of  statistics  we  know  that  the  optimal 
weighting  matrix  for  this  purpose  is  the  inverse  of  the  data 
covariance  matrix.  Therefore  instead  of  minimizing  (41),  we 
minimize 

€7C'‘§  =f^2-^)TC'1(A?-'V) 
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where  C  is  the  covariance  matrix  of  data.  Differentiating 
(42)  with  respect  to  XT  and  setting  it  equal  to  zero  gives 

x  =  (atc'1a)‘1x,c15 


Thus  we  have  obtained  the  best  estimate  X  of  the  true  solu¬ 
tion  X  for  a  given  level  of  error  in  the  data. 

Hamilton  (1964)  has  shown  that  under  the  above  circum¬ 
stances  the  following  holds: 

covcj)  =  CATAy‘  q1  m) 


where  cov(X)  is  the  covariance  matrix  of  the  model  para¬ 
meters,  and  is  the  variance  of  the  data.  Substituting 

(25)  into  (44)  for  our  generalized  linear  inverse  formalism 
gives 

C0V(g).  Vt  A'j  w  (4S) 


Hence  the  variance  of  our  model  parameters  may  be  expressed 
as 


(46) 


8.  Consideration  for  Model  Parameters  with  Different  Dimensions 
Since  in  the  present  problem,  we  have  model  parameters 
pC  ,  |3  .  f3  and  oi  which  have  different  units  and  sizes, 
we  wish  to  account  for  this  property  of  our  set  of  model 
parameters . 
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Smith  and  Franklin  (1969)  have  shown  that  the  generalized 
linear  inversion  of  the  system  (10)  simultaneously  minimizes 


both 


6 T  €  =  A*-H)  -|  Ax- v 


(41) 


XT  x  =  ||Xjj2 


(41) 


In  other  words ,  the  generalized  linear  inverse  selects  the 
smallest  variation  X  of  the  parameters  that  will  satisfy 
the  system  of  simultaneous  linear  equations  (10) . 

In  order  to  account  for  the  different  dimensions  or  sizes 
of  the  model  parameters ,  we  introduce  the  inverse  of  the 
parameter  covariance  matrix,  ' ,  as  a  weighting  matrix. 

Wiggins  (1972)  has  shown  that  this  choice  of  weighting  has 
the  interesting  and  convenient  property  of  making  the  lengths 
of  the  rows  of  AW  invariant  to  the  selection  of  layer  thick¬ 
nesses.  Now  we  want  to  simultaneously  minimize 


€  C  € 


(4D 


x  W'Ax 


(48) 


instead  of  (41)  and  (47) . 

In  order  to  simplify  the  process,  we  introduce  the 
transformation  of  Wiggins  (1972)  and  Jackson  (1972) : 

4*sC‘V:A  Wi/2 

x*  =  yi/2  x 


(49) 
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y  *  C'"1  H  <«> 

When  we  consider  the  system 

A*  x*  =  (Si.) 

£’  =  A*  2*  -  3“  U3) 

we  find  that 

6*T  6*  =  €t  C  1  §  (54) 

2*Tx*  =  xT  W  '*1  f55) 


Thus  we  have  shown  that  with  transformations  (49)  to  (51) , 
minimizing  £*  €*  and  £  &  is  equivalent  to  minimizing 

€T£_1€  and  XTh/"*X-  Therefore  we  take  the  convenience  of 
doing  the  former. 


RESULTS  AND  DISCUSSION 


The  corrections  for  the  phase  response  of  both  the  old 
and  new  instruments  are  shown  in  Table  3.  Both  the  period 
and  corrections  are  in  units  of  seconds. 

Der  et  al.  (1970)  have  shown  that  when  group-velocity 
dispersion  observation  points  are  evenly  spaced  on  a  logarith¬ 
mic  period  scale,  each  of  the  chosen  data  points  furnishes 
approximately  the  same  amount  of  information  about  the  sub¬ 
surface.  In  other  words,  this  choice  of  observation  periods 
provides  the  greatest  independence  of  data  with  the  same  num¬ 
ber  of  observation  points.  In  the  present  study,  observation 
points  were  chosen  at  approximately  equally  spaced  intervals 
on  a  logarithmic  period  scale  with  the  interval  such  that 
A  (log  T)  <=55  0.1  in  the  range  from  10  to  100  sec.  This 
choice  of  interval  or  the  number  of  observation  points  was 
guided  by  inspection  of  the  shape  of  partial-derivative  curves 
obtained  in  preliminary  studies  of  the  structure  of  the  area. 
The  periods  indicated  in  Table  3  represent  those  observation 
periods  so  determined. 

The  corrections  for  the  sphericity  and  gravity  of  the 
earth  are  shown  in  Table  4.  The  values  in  Taole  4  were 
weighted  in  proportion  to  the  path  lengths  corresponding 
to  different  provinces  as  indicated  in  Table  2;  they  were  sub¬ 
tracted  from  the  observed  group  velocities.  As  noted  from 
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Table  4,  Che  corrections  are  so  small  that  they  are  less  than 
the  uncertainty  of  the  observed  data  discussed  previously. 
However,  since  they  provide  a  systematic  change  in  group 
velocity,  these  corrections  were  applied  to  the  observed  data. 

The  observed  group  velocity  dispersion  curves  determined 
from  the  phase-matched  filtering  and  corrected  as  discussed 
above  are  shown  in  Fig.  5.  In  comparison  with  dispersion 
curves  of  an  average  earth  (e.g.  Oliver,  1962;  to  avoid  con¬ 
fusion  due  to  crowded  curves,  they  are  not  shown  in  Fig.  5),  the 
following  is  observed.  The  dispersion  curve  obtained  for 
continental  Alaska  (event  0509)  is  similar  to  that  of  the 
average  continental  earth  of  Oliver  (1962)  except  that  it  has 
slightly  higher  velocities  in  the  longer  periods  and  lower 
velocities  in  the  shorter  periods.  For  epicenters  farther 
to  the  west  along  the  Aleutian  Islands,  the  dispersion  curves 
shift  toward  that  of  the  average  oceanic  earth  of  Oliver 
(1962)  but  the  velocity  values  are  still  much  closer  to  that 
of  the  average  continental  earth  than  to  the  average  oceanic 
earth . 

The  results  of  the  inversion  of  the  group-velocity  dis¬ 
persion  relationships  are  described  in  the  following 
paragraphs.  In  the  inversion  scheme  used  here,  since  there 
is  no  way  to  predict  precisely  the  variance  of  parameters, 
a  constant  value  of  0.15  was  used  as  an  initial  estimate  of 
the  standard  deviation  of  parameters. 
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1 .  Continental  Alaska 

The  initial  model  for  the  structure  of  continental 
Alaska  is  shown  in  Table  5.  In  constructing  initial  models, 
the  updated  comprehensive  parametric  earth  models  of 
Dziewonski  e£  al.  (1975)  were  frequently  consulted.  Stacey 
(1977)  gives  these  models  as  representative  earth  models  in 
his  recent  textbook. 

Although  the  method  developed  in  this  study  and  the 
computer  code  used  to  implement  the  method  are  capable  of  com¬ 
puting  the  group-velocity  partial  derivatives  with  respect  to 
thickness  of  the  layers,  for  the  sake  of  simplicity  and  con¬ 
venience,  the  thicknesses  of  the  layers  were  held  fixed  and 
parameters  fi  ■  f  and  oi  within  each  layer  were  corrected 
in  the  inversion  process. 

Partial  derivatives  of  the  group  velocity  with  respect 
to  ,  p  and  oC  are  shown  in  Figs.  6,  7  and  8,  respectively. 
The  numbers  on  curves  designate  the  corresponding  layers. 

It  is  noted  from  the  figures,  as  expected,  that  partial  deri¬ 
vative  maxima  shift  toward  longer  periods  as  the  depths  of 
the  layer  increase.  It  is  also  noted  that  the  amplitudes  of 
the  ft  partial  derivatives  are  significantly  greater  than 
those  of  the  other  partials  and  that  the  oC  partials  are  the 
least  significant  of  all. 

The  structure  of  the  final  model  for  this  region  obtained 
from  the  Alaskan  event  is  shown  in  Fig.  9.  Also  shown  in  the 
figure  are  standard  deviations  of  the  model  parameters 
indicated  by  horizontal  bars  at  the  center  of  each  layer. 
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TABLE  5 


INITIAL 

MODEL  FOR  CONTINENTAL 

ALASKA 

Layer 

P-velocity 

(km/sec) 

S-velocity 
(km/ sec) 

Density 

(g/cm3) 

Thickness 

(km) 

Depth  to 
Bottom  of 
Layer  (km) 

1 

3.670 

2.310 

2.320 

0.8 

0.8 

2 

5.411 

3.266 

2.637 

2.1 

2.9 

3 

5.789 

3.453 

2.768 

6.1 

9.0 

4 

6.152 

3.483 

2.800 

7.0 

16.0 

5 

6.398 

3.617 

2.848 

7.0 

23.0 

6 

6.490 

3.894 

2.962 

10.0 

33.0 

7 

6.687 

4.024 

2.984 

10.0 

43.0 

8 

8.016 

4.649 

3.325 

20.0 

63.0 

9 

8.021 

4 . 666 

3.333 

50.0 

113.0 

10 

7.853 

4.461 

3.387 

50.0 

163.0 

11 

7.853 

4.468 

3.388 

50.0 

213.0 

12 

8.651 

4.652 

3.452 

Figure  6:  Group  Velocity  Partial  Derivatives  with  Respect 


Figure  7:  Group  Velocity  Partial  Derivatives  with  Respect 
to  Density  in  Alaska 


DEPTH  (KM) 
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Figure  9:  Structure  of  the  Final  Model  for 
Continental  Alaska 
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The  greatest  estimated  standard  deviation  of  parameters 
was  0.28  km/ sec  for  the  shear  velocity  of  the  layer  centered 
at  a  depth  of  188  km.  The  observed  and  predicted  group  velo¬ 
cities  are  shown  in  Fig.  10.  The  rms  error  of  the  predicted 
values  against  the  observed  ones  was  0.02  km/sec,  which  is  less 
than  the  uncertainty  of  the  data.  The  differences  between  the 
observed  and  predicted  values  were  all  within  the  standard 
deviation  of  the  data  except  at  the  period  of  25  sec,  where 
the  difference  was  0.04  km/ sec. 

It  is  observed  from  Fig.  9,  and  also  from  Figs.  16  and 
23  for  the  other  regions  as  shown  later,  that  the  magnitude 
of  the  standard  deviations  of  model  parameters  decreases  from 
j3  to  P  and  to  oL  .  These  trends  do  not  appear  to  represent 
the  real  accuracy  of  the  model,  for  from  the  physical  point  of 
view,  would  have  about  the  same  as  or  possibly  greater  un¬ 
certainty  than  /3  .  Therefore  this  apparent  differential 
accuracy  among  the  parameters  is  interpreted  to  be  caused  by 
some  property  of  the  present  generalized  linear  inversion 
method.  A  possible  explanation  is  given  below.  As  an  initial 
estimate  of  standard  deviation  of  parameters,  we  used  0.15 
for  all  parameters.  Since  the  values  of  partial  derivatives 
of  the  group  velocity  with  respect  to  oC  are  much  smaller  than 
those  of  ,  and  still  smaller  than  those  of  /5  ,  while  a 
constant  initial  estimate  of  standard  deviation  was  used  for 
all  parameters,  less  flexibility  in  parameter  correction  seems 
to  have  been  given  to  f*  and  still  less  to  oC  in  order  for  the 
system  to  find  a  new  model  which  fits  the  observed  data  better 
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Figure  10:  Observed  and  Predicted 
Group  Velocities  for 
Continental  Alaska 


than  the  initial  model.  In  other  words,  smaller  allowance 
or  variance  for  j*  and  relative  to  was  required  by  the 
system  in  making  correction  to  the  model.  Therefore  the 
apparent  smaller  standard  deviations  for  and  oC  are  attri¬ 
buted  to  estimation  error  rather  than  to  greater  accuracy  in 
determining  parameters  />  and  oi  •  A  better  initial  estimate 
of  the  covariance  matrix  of  model  parameters  may  resolve  this 
problem. 

Two  features  of  the  structure  worth  noticing  from  Fig.  9 
are  the  depth  to  the  Mohorovicic  discontinuity  and  the  pre¬ 
sence  of  the  low  velocity  zone  (LVZ) .  The  depth  to  the  Moho 
is  43  km.  Although  no  previous  studies,  to  this  author's 
knowledge,  have  been  conducted  covering  the  same  area  as  the 
present  one,  the  thickness  of  the  crust  obtained  here  may  be 
compared  with  results  of  other  studies  made  in  areas  nearby. 
Using  gravity  data,  Woollard  et  al.  (1960)  estimated  the 
crustal  thickness  near  Fairbanks  to  be  about  33  km.  Hales  and 
Asada  (1966)  gave  48  to  53  km,  estimated  from  their  seismic 
refraction  studies,  as  the  depth  to  the  Moho  for  the  area 
from  College  Fiord  to  the  northeast,  which  is  near  Fairbanks. 
Hanson  et  al.  (1968)  found  from  their  seismic  refraction 
interpretation  that  the  Moho  is  at  a  depth  of  32  km  under 
Fairbanks  and  48  km  under  their  shot  point,  which  is  about 
120  km  southwest  of  Fairbanks.  The  present  model  is  essen¬ 
tially  in  agreement  with  the  refraction  results  of  Hanson 
et  al.  (1968)  in  the  structure  of  the  top  10  km,  except  that 
their  model  has  a  layer  of  sediments  a  little  thicker  than  in 
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the  present  model.  Since  their  model  has  only  a  few  layers, 
the  present  model  cannot  be  checked  with  their  results  in 
the  deeper  structure. 

From  Fig.  9  it  is  also  noted  that  a  LVZ  is  present  in 
the  region  and  it  extends  from  the  depth  of  about  113  to  213 
km.  The  presence  and  depth  range  of  the  LVZ  of  the  present 
model  is  in  essential  agreement  with  the  parametric  conti¬ 
nental  model  of  Dziewonski  et  al.  (1975) . 

Resolving  kernels  for  ft  and  of  the  final  model  for 

this  region  are  shown  in  Figs.  11  and  12,  respectively.  Depth 
corresponding  to  each  resolving  kernel  is  indicated  in  the 
figures.  Some  sidelobes  and  deviation  of  the  resolution 
matrix  from  an  identity  matrix  are  observed  in  the  resolving 
kernels  for  ft  .  The  sidelobes  and  deviation  are  seen  to  have 
worsened  in  the  case  of  P  .  The  resolving  kernels  for  oC  in 
this  and  the  other  regions  were  so  irregular  in  shape  that 
they  were  not  plotted. 

Also  shown  for  the  bottom  layer  of  the  crust  in  Fig.  11 
is  the  spread,  which  is  the  depth  range  of  the  rectangle  which 
has  the  same  area  as  that  under  the  associated  resolving  ker¬ 
nel.  Since  the  spread  for  this  layer  is  16  km  compared  with 
a  10  km  thickness  of  the  layer,  the  uncertainty  of  the  depth 
to  the  Moho  is  interpreted  to  be  +  3  km.  Therefore  a  reason¬ 
able  value  for  the  thickness  of  the  crust  in  this  region  is 
43  +  3  km,  where  43  km  is  the  depth  to  the  bottom  of  the 
lowest  layer  in  the  crust.  Although  the  corresponding  spread 
for  P  is  comparable  to  that  for  {3  ,  since  (3  is  the  most 
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Figure  11:  Resolving  Kernels  for  Shear  Velocity 
in  Continental  Alaska 
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Figure  12:  Resolving  Kernels  for  Density 
in  Continental  Alaska 
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influential  parameter  in  the  inversion  and  its  resolving  ker¬ 
nels  look  most  reliable  in  the  sense  that  they  show  least 
sidelobes,  the  spread  is  shown  only  for  (3  . 

The  information  density  distributed  in  the  model  for  this 
region  is  shown  in  Fig.  13.  It  is  noted  from  the  figure  that 
the  information  density  matrix  has  to  some  extent  been  degra¬ 
ded,  i.e.  distorted  from  the  identity  matrix.  This  degree  of 
distortion  is,  however,  considered  desirable.  If  the  informa¬ 
tion  density  matrix  is  a  perfect  identity  matrix,  although  the 
economy  of  data  acquisition  and  processing  is  optimal,  there 
always  exists  a  possibility  of  missed  information  needed  in 
the  inversion.  On  the  other  hand,  if  the  information  density 
matrix  is  too  severely  distorted,  it  may  guarantee  that  suffi¬ 
cient  information  may  have  been  propagated  into  the  inversion 
process,  but  the  economy  of  data  acquisition  and  processing 
suffers  severely.  Although  the  desirable  degree  of  distortion 
of  the  information  density  matrix  is  difficult  to  determine 
quantitatively  and  consequently  is  a  subjective  matter,  that 
degree  of  distortion  which  appears  in  Fig.  13  is  considered  a 
reasonable  compromise. 

2.  Bering  Shelf 

The  initial  model  for  the  Bering  Shelf  region  is  shown 
in  Table  6.  Group-velocity  partial  derivatives  with  respect 
to  the  three  parameters  for  this  region  are  shown  for  periods 
25  and  40  sec  in  Figs.  14  and  15,  respectively.  As  expected, 
the  partial  derivative  maxima  and  minima  are  shifted  to 
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TABLE  6 

INITIAL  MODEL  FOR  THE  BERING  SHELF 


.aver 

HuSSSSl 

y  Density 
(g/cm3) 

Thickness 

(km) 

Depth 
;  Bottoi 
Layer 

0 

1 

2 

3 

1.520 

2.000 

3.670 

5.789 

0.000 

1.000 

2.310 

3.453 

1.030 

1.500 

2.320 

2.768 

0.1 

1.0 

2.3 

4.2 

0.1 

1.1 

3.4 

7.6 

4 

6.152 

3.483 

2.800 

7.0 

14.6 

5 

6.393 

3.617 

2.848 

7.0 

21.6 

6 

6.490 

3.894 

2.962 

6.0 

27.6 

7 

7.934 

4.654 

3.310 

10.0 

37.6 

8 

8.016 

4.649 

3.325 

20.0 

57.6 

9 

8.021 

4.666 

3.333 

50.0 

107.6 

10 

7.853 

4.461 

3.387 

50.0 

157.6 

11 

7.853 

4.468 

3.388 

55.0 
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8.651 

4.652 
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PARTIAL  DERIVATIVES 


Figure  14:  Group  Velocity  Partial  Derivatives 
for  1=25  in  the  Bering  Shelf 


greater  depths  as  the  period  increases.  It  is  also  noted 
from  these  figures  and  Figs.  21  and  22  as  shown  later,  that 
the  amplitude  of  density  partial  derivatives  relative  to  that 
of  shear  velocity  partial  derivatives  decreases  with  depth. 

The  structure  of  the  final  model  for  this  region  derived 
from  the  analysis  of  event  1410,  along  with  the  standard 
deviations  of  the  parameters,  is  shown  in  Fig.  16.  The 
greatest  standard  deviation  was  0.29  km/sec  for  shear  velocity 
of  the  layer  centered  at  a  depth  of  133  km.  The  observed  and 
predicted  group  velocities  are  shown  in  Fig.  17.  All  the 
differences  between  the  observed  and  predicted  values  and 
their  rms  error  were  within  the  uncertainty  of  the  data. 

It  is  seen  from  Fig.  16  that  the  depth  to  the  Moho  is 
27.6  km.  This  value  may  be  compared  with  29  km  estimated 
from  seismic  refraction  studies  of  Shor  (1964)  in  the  south¬ 
western  part  of  the  Bering  Shelf  (stations  MK11  and  MK12  of 
Shor,  1964).  Also  the  thicknesses  and  velocities  of  surface 
layers  are  in  reasonable  agreement  with  those  of  Shor  (1964) . 
Since  Shor's  model  has  only  a  few  layers,  the  deeper  layers 
cannot  be  checked  against  the  refraction  results.  Also  noted 
from  Fig.  16  is  the  existence  of  a  LVZ  which  ranges  from  about 
108  km  to  about  213  km  in  depth.  This  depth  range  of  the  LVZ 
may  be  comparable  to  the  earth  model  averaged  over  the  whole 
earth  including  oceans  and  continents,  of  Dziewonski  et  al . 
(1975). 

Resolving  kernels  for  {$  and  along  with  the  corres¬ 
ponding  depths  of  the  present  model  for  this  region  are  shown 
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in  Figs.  18  and  19,  respectively.  The  weaker  resolving 
power  of  p  ,  and  the  still  weaker  resolving  power  of  oC 
which  is  not  shown  graphically,  compared  with  that  of  @ 
are  noted  from  these  figures  as  well  as  from  Figs.  11,  12, 

25  and  26.  In  Fig.  18,  the  spread  is  shown  for  the  bottom 
layer  of  the  crust.  The  spread  is  seen  to  be  14  km.  A 
reasonable  estimate  for  the  thickness  of  the  crust  for  this 
region  is  28  +  4  km. 

The  information  density  of  the  group  velocities  used  in 
the  inversion  process  for  the  Bering  Shelf  region  is  shown  in 
Fig.  20.  It  is  seen  to  have  been  distorted  from  the  identity 
matrix  to  about  the  same  extent  as  in  the  continental  Alaska 
region . 

3.  Aleutian  Basin 

The  initial  model  for  the  Aleutian  Basin  region  is  shown 
in  Table  7.  The  group-velocity  partial  derivatives  with  res¬ 
pect  to  the  earth  parameters  for  this  region  are  shown  for 
periods  16  and  50  sec  in  Figs.  21  and  22,  respectively.  Since 
the  two  periods  represented  in  these  figures  are  farther  apart 
than  in  Figs.  14  and  15  for  the  Bering  Shelf  region,  the  depth 
shift  of  the  partial-derivative  maxima  is  seen  to  be  more  pro¬ 
nounced  here  than  in  Figs.  14  and  15.  It  is  also  clearly  seen 
from  these  figures  that  as  the  period  of  waves  increases,  the 
shape  of  the  partial  derivatives  broadens  accordingly. 

The  final  model  obtained  by  the  analysis  of  event  2224, 
with  standard  deviations  of  the  parameters,  is  shown  in  Fig.  23. 
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Figure  18:  Resolving  Kernels  for  Shear 
Velocity  in  the  Bering  Shelf 
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TABLE  7 


INITIAL  MODEL  FOR  THE 

ALEUTIAN 

BASIN 

P-velocity 

S-velocity 

Density 

Thickness 

Depth  to 
Bottom  of 

Layer 

(km/sec) 

(km/sec) 

(ft/cm3} 

(km) 

Layer  (km) 

0 

1.520 

0.000 

1.030 

3.0 

3.0 

1 

2.000 

1.000 

1.500 

1.0 

4.0 

2 

4.556 

2.500 

2.480 

4.8 

8.8 

3 

6.333 

3.550 

2.867 

2.5 

11.3 

4 

7.000 

4.000 

3.000 

6.5 

17.8 

5 

7.900 

4.550 

3.305 

40.0 

57.8 

6 

7.873 

4.335 

3.359 

40.0 

97.8 

7 

7.873 

4.335 

3.383 

40.0 

137.8 

8 

7.873 

4.335 

3.426 

70.0 

207.8 

9 

7.873 

4.335 

3.432 

10.0 

217.8 

10 

8.651 

4.652 

3.452 

Figure  21:  Group  Velocity  Partial  Derivatives 
for  T=16  in  the  Aleutian  Basin 
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Figure  23:  Structure  of  the  Final  Model  for  the 
Aleutian  Basin 
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Fig.  24.  All  the  differences  between  the  observed  and  pre¬ 
dicted  values  and  their  rms  error  were  within  the  uncertainty 
of  the  data.  It  is  found  from  Fig.  23  that  the  crustal  thick¬ 
ness  is  17.8  km.  Shor  (1964)  has  estimated  the  depth  to  the 
Moho  in  this  region  as  14  km,  which  was  obtained  from  his 


refraction  analysis.  From  the  results  of  studies  on  travel 
times  and  waveforms  of  both  refracted  and  reflected  body  waves, 
Helmberger  (1968)  has  presented  models  of  the  Aleutian  Basin. 
His  model  with  no  transition  zone  between  the  crust  and  man¬ 


tle  gives  19  km  (his  Fig.  10  for  station  L13)  as  the  depth  to 
the  boundary,  while  his  model  with  one  transition  layer  of 
thickness  less  than  1  km  shows  the  depth  of  the  crust- to-mantle 
transition  layer  around  16  to  17  km  (his  Fig.  12). 

While  Helmberger' s  model  with  the  thick  transition  zone 
for  his  station  L9  may  be  applicable  to  that  part  of  the 
Aleutian  Basin  which  is  near  the  Aleutian  Islands,  the  present 
result  of  about  18  km  may  represent  the  average  structure  of 


the  Aleutian  Basin.  It  may  well  be  pointed  out  that  the 
great  circle  path  of  event  2224  traverses  the  middle  of  the 
Basin.  The  crustal  structure  of  the  present  model  is  in 
reasonable  agreement  with  nelmberger's  model  (his  Fig.  10). 

The  present  model  in  Fig.  23  also  shows  the  presence  of 
a  LVZ  which  extends  from  the  depth  of  about  60  km  to  about 
220  km.  The  presence  and  depth  range  of  the  present  model 
agree  with  the  parametric  oceanic  earth  model  of  Dziewonski 
et  al.  (1975) . 
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The  resolving  kernels  of  yS  and  of  the  final  model 
for  this  region  are  shown  in  Figs.  25  and  26,  respectively. 
Depths  associated  with  each  kernel  are  also  shown  in  the  fig¬ 
ures.  The  spread  for  the  deepest  crustal  layer  is  also  shown 
in  connection  with  the  resolving  kernel  for  ft  of  the  layer. 

The  spread  is  14  km.  A  reasonable  estimate  for  the  crustal 
thickness  of  the  Aleutian  Basin  is  18+4  km.  Shor's  value 
of  14  km,  which  is  smaller  than  the  present  and  Helmberger's 
estimates,  may  be  regarded  as  a  lower  limit  of  the  depth  to 
the  Moho . 

The  present  result  on  the  structure  of  the  Aleutian  Basin 
shows  that,  although  its  crustal  thickness  is  greater  than  in 
average  oceans,  this  part  of  the  Bering  Sea  is  oceanic  and 
has  not  been  continentalized  in  the  characteristics  of  the 
crust.  However,  its  possibility  of  conversion  from  an  ocean 
basin  into  a  continental  mass  in  the  future  as  noted  by  Shor 
(1964)  remains  to  be  studied. 

The  information  density  of  the  observed  group-velocity 
dispersion  relationship  used  in  the  inversion  for  this  region 
is  shown  in  Fig.  27.  It  is  noted  from  the  figure  that  the 
information  density  matrix  for  this  case  has  been  degraded 
to  a  greater  extent  than  in  cases  of  continental  Alaska  and 
the  Bering  Shelf.  The  cause  of  this  greater  distortion  is 
thought  to  be  due  to  the  simpler  shape  of  the  dispersion  curve 
for  this  region  compared  with  those  of  the  other  regions, 
while  the  number  of  observation  periods  was  the  same  in  all 


cases  . 
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Figure  25:  Resolving  Kernels  for  Shear  Velocity  in 
the  Aleutian  Basin 
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Figure  26:  Resolving  Kernels  for  Density  in  the  Aleutian 
3asin 


Figure  27:  Information  Density  in  the  Aleutian  Basin 
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Thus  we  have  obtained  estimates  of  the  structure  of  the 
three  provinces  of  the  study  area.  In  order  to  utilize  the 
additional  data  available,  the  observations  from  the  events 
1519  and  2005,  and  to  check  the  results  obtained  above,  the 
following  efforts  were  made.  From  the  experience  gained  to 
this  point,  the  author  has  learned  that  it  took  many  trials 
and  accompanying  time  of  considerable  amount  to  find  a  "good" 
initial  model  which  gives  a  satisfactory  convergence  of  pre¬ 
dicted  data  to  the  observed  one.  Therefore  instead  of 
repeating  the  inversion  procedures  with  the  rest  of  the 
observed  data,  it  was  decided  to  solve  a  forward  problem  with 
the  structures  obtained  and  with  appropriate  proportions  of 
path  length  and  to  compare  the  results  with  observed  data. 

With  the  information  as  shown  in  Figs.  9,  16  and  23, 
and  Table  2,  group  velocities  were  calculated  for  events 
1519  and  2005  and  the  results  along  with  the  observed  group 
velocities  are  shown  in  Figs.  28  and  29.  As  can  be  seen  from 
the  figures,  the  agreement  between  the  observed  and  predicted 
values  was  very  good  except  at  a  few  of  the  shortest  periods. 
These  results  are  therefore  in  agreement  with  the  final 
models  we  have  obtained. 
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Figure  29: 
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CONCLUSION 


From  the  studies  of  surface  waves  by  the  method  des¬ 
cribed  in  this  study  and  the  results  discussed  in  the  last 
chapter,  the  following  conclusions  are  drawn: 

1.  Phase-matched  filtering  has  been  effective  in 
determining  group-velocity  dispersion  curves  in  the  present 
study  area.  The  curves  have  been  inverted  and  checked  with 
previous  results  of  gravitational,  seismic  refraction  and 
reflection  studies. 

2.  The  best  estimates  of  the  group-velocity  dispersion 
relationships  in  the  area  studied  have  been  obtained  and 
presented  in  Fig.  5. 

3.  The  great  circle  paths  for  the  events  used  in  this 
study  appear  to  coincide  with  the  least-time  paths.  This  was 
the  case  due  to  the  geometric  arrangements  of  the  physio¬ 
graphic  provinces  and  their  boundaries  of  the  study  area 
relative  to  the  source-station  pairs. 

4.  The  method  of  studying  a  compound  area  by  dividing 
it  into  two  or  more  "pure-path"  provinces  proved  effective 
under  the  circumstances  described  in  the  last  sentence.  The 
group-velocity  dispersion  relationships  derived  from  this 
method  for  the  "pure-path"  provinces  of  the  Bering  Shelf 
and  the  Aleutian  Basin  are  presented  in  Figs.  17  and  24, 
respectively . 
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5.  The  exact  analytical  method  formulated  ir  this  study 
for  the  computation  of  Rayleigh-wave  phase-velocity  partial 
derivatives  with  respect  to  earth  parameters  has  proved 
effective,  has  not  shown  any  numerical  difficulties  when 
coded  with  double  precision,  and  has  been  successfully  incor¬ 
porated  into  a  generalized  linear  inversion  algorithm. 

6.  A  generalized  linear  inversion  method  has  been  applied 
to  an  inversion  of  group-velocity  dispersion  relationships  in 
order  to  find  the  shear-wave  velocity,  density  and  compressional- 
wave  velocity  structure  of  the  earth.  The  results  have  yielded 
the  following  estimates: 

(i)  The  average  crustal  thickness  of  the  continen¬ 
tal  Alaska  region  studied  is  43  +  3  km.  A  LVZ  is  present  in 
this  region  and  extends  from  about  113  to  about  213  km  depth. 

(ii)  The  depth  to  the  Mono  in  the  Bering  Shelf 
region  is  28  +  4  km.  A  LVZ  ranges  in  depth  from  about  108 
to  about  213  km. 

(iii)  The  crust  in  the  Aleutian  Basin  region  is 
oceanic  in  nature  and  the  boundary  between  the  continental 
Bering  Shelf  and  the  oceanic  Aleutian  Basin  appears  to  lie 
near  and  parallel  to  the  1000-fathom  bathymetric  contour  line. 

(iv)  The  thickness  of  the  crust  in  the  Aleutian 
Basin  region  is  18+4  km.  A  LVZ  is  present  in  this  region 
and  extends  from  about  60  to  about  220  km  depth. 


APPENDIX 


Expressions  necessary  to  compute  the  phase-velocity  par¬ 
tial  derivatives  are  given  below.  It  is  hoped  that  these 
expressions,  although  somewhat  lengthy,  will  give  those  who 
use  the  method  of  this  paper  a  convenience  and  saving  of  time 
required  to  carry  through  the  very  time-consuming  and  tedious 
calculations.  It  should  be  noted  that  in  the  case  of  Rayleigh 
waves,  there  are  three  situations  in  connection  with  the  rela¬ 
tion  between  the  phase  velocity  and  body-wave  velocities . 

They  are:  (1)  c>o(m  and  c>^m,  i.e.,  both  m  and  Jl ^  m  are 
positive  real;  (2)  c<o^m  and  c>  P  m  >  1  • e  •  >  -/Lot.  m  is  negative  pure 
imaginary  while  is  positive  real;  (3)  c<p(m  and  c<  P  m> 

i.e.,  both  and  /L £m  are  negative  pure  imaginary.  In  the 

following  expressions  a  triple  sign--three  signs  put  together 
vertically--applies ,  in  the  order  from  top  to  bottom,  to  the 
above  three  cases.  When  the  three  cases  have  a  common  sign, 
only  a  single  sign  appears.  It  is  also  understood  that  when¬ 
ever  and  are  pure  imaginary,  and  in  the 

following  expressions  actually  represent  m  an<^  m  , 

respectively;  and  trigonometric  functions  become  corresponding 
hyperbolic  functions. 

Since  the  elements  of  the  layer  matrix  have  been  published 
previously,  see  Haskell  (,1953,  p.  21)  and  Harkrider  (1964,  Eq . 
16),  only  their  partial  derivatives  are  presented  here. 
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